The present invention relates to the field of signal processing, and in particular, to an improved system and method for computing the amount of shift between two signals.
A great variety of systems may be characterized by the property that they generate a series of signals which are shifted and noise-perturbed versions of some known template signal. For example, radar systems operate by repeatedly transmitting a known pulse template and receiving the pulse reflections generated by objects (such as aircraft) which intercept the field of view of the radar beam. The pulse reflections are delayed versions of the transmitted pulse and are corrupted by various forms of noise such as multipath noise. The time-delay between the transmitted pulse and the pulse reflection is a measure of an object""s distance.
The problem of estimating the shift (e.g. time-delay) between two signals x(n) and y(n) is one of the most fundamental problems in signal processing. Perhaps the most popular method for solving the shift estimation problem involves the use of cross-correlation. The cross-correlation method involves computing a cross-correlation                                                         R              xy                        ⁡                          (              τ              )                                =                                    ∑              n                        ⁢                                          x                ⁡                                  (                                      n                    +                    τ                                    )                                            ⁢                              y                ⁡                                  (                  n                  )                                                                    ,                            (A1)            
for shift values xcfx84 spanning a range of expected values. The shift-value xcfx84max for which the cross-correlation Rxy (xcfx84) is maximized is taken to be an estimate for the shift between signals x(n) and y(n).
The problem of shift estimation is just as relevant for periodic signals, or more generally, for signals defined on the integers modulo N, where N is a positive integer. (Physically realizable discrete-time signals are finite in duration and may be naturally interpreted as being defined on the integers modulo N.) Given two signals f and g defined on the integers modulo N, where signal g is known to be a shifted version of signal f the circular cross-correlation                                           R            fg                    ⁡                      (            d            )                          =                              ∑                          n              =              0                                      N              -              1                                ⁢                                    f              ⁡                              (                                                      [                                          n                      +                      d                                        ]                                    N                                )                                      ⁢                          g              ⁡                              (                                                      [                    n                    ]                                    N                                )                                                                        (A2)            
may be evaluated for shift value d equal to 0, 1, 2, . . . , Nxe2x88x921. The shift-value dmax for which the circular cross-correlation Rfg (d) is maximized may be taken as an estimate for the circular shift between signals f and g. Thus, the circular cross-correlation method requires N2 multiplies and N(Nxe2x88x921) additions in order to generate a shift estimate between two functions. More efficient methods do exist. All of them are complicated and cannot reduce the numerical complexity significantly.
It can be shown that the circular cross-correlation of signals f and g may be expressed as a circular convolution, i.e.                                                         R              fg                        ⁡                          (              d              )                                =                                    ∑                              k                =                0                                            N                -                1                                      ⁢                                          f                ⁡                                  (                                                            [                                              d                        -                        k                                            ]                                        N                                    )                                            ⁢                              u                ⁡                                  (                                                            [                      k                      ]                                        N                                    )                                                                    ,                            (A3)            
where signal u denotes the time-reversal of signal g, i.e.
u([k]N)=g([xe2x88x92k]N).xe2x80x83xe2x80x83(A4)
Thus, the discrete Fourier transform {circumflex over (R)}fg of circular cross-correlation function Rfg equals the product of the discrete Fourier transforms {circumflex over (f)} and û of signals f and u respectively, i.e.
xe2x80x83{circumflex over (R)}fg={circumflex over (f)}xc2x7û.xe2x80x83xe2x80x83(A5)
Thus, the circular cross-correlation function Rfg may be computed indirectly as suggested by the following expression:
Rfg=IDFT({circumflex over (f)}xc2x7ĝ),xe2x80x83xe2x80x83(A6)
where IDFT denotes the inverse discrete Fourier transform operator. If Fast Fourier Transforms (FFTs) are used to compute the pair of forward DFTs and the inverse DFT, this method for computing the circular cross-correlation function Rfg has computational complexity of order xcex1N log2(N), where xcex1 is a fixed number. Thus, this indirect method for computing the circular cross-correlation function is significantly more efficient than the exhaustive cross-correlation method as represented by equation (A2). However, according to Mitra et al., it is only useful for N greater than 256, due to complex number overhead and due to the factor xcex1.
Given a signal h defined on the integers modulo N, the Discrete Fourier transform ĥ of signal h may be defined as                                                         h              ^                        ⁡                          (              k              )                                =                                    ∑                              n                =                0                                            N                -                1                                      ⁢                                          h                ⁡                                  (                  n                  )                                            ·                              exp                ⁡                                  (                                                                                    -                        j2                                            ⁢                                              xe2x80x83                                            ⁢                      π                      ⁢                                              xe2x80x83                                            ⁢                      nk                                        N                                    )                                                                    ,                            (A7)            
where index k varies from 0 to Nxe2x88x921. Let fd denote the signal obtained by shifting the elements of signal f by d positions, i.e. fd([n]N)=f([n+d]N). Using the DFT definition (A7), it can be shown that the discrete Fourier transform {circumflex over (f)}d of shifted signal fd is given by                                                                         f                ^                            d                        ⁡                          (              k              )                                =                                                    f                ^                            ⁡                              (                k                )                                      ·                          exp              ⁡                              (                                                      j                    ⁢                                          xe2x80x83                                        ⁢                    dk2                    ⁢                                          xe2x80x83                                        ⁢                    π                                    N                                )                                                    ,                            (A8)            
where k varies from 0 to Nxe2x88x921, and where {circumflex over (f)} denotes the discrete Fourier transform of signal f Expression (A8) implies that                                                         Arg              ⁡                              [                                                                            f                      ^                                        d                                    ⁡                                      (                    k                    )                                                  ]                                      -                          Arg              ⁡                              [                                                      f                    ^                                    ⁡                                      (                    k                    )                                                  ]                                              =                                                    dk2                ⁢                                  xe2x80x83                                ⁢                π                            N                        +                          2              ⁢                              xe2x80x83                            ⁢              π              ⁢                              xe2x80x83                            ⁢              r                                      ,                            (A9)            
for some integer r, where Arg[z] denotes the principle angle of complex number z. Solving for shift parameter d gives the expression                     d        =                                            N                              2                ⁢                                  xe2x80x83                                ⁢                π                ⁢                                  xe2x80x83                                ⁢                k                                      ⁢                          {                                                Arg                  ⁡                                      [                                                                                            f                          ^                                                d                                            ⁡                                              (                        k                        )                                                              ]                                                  -                                  Arg                  ⁡                                      [                                                                  f                        ^                                            ⁡                                              (                        k                        )                                                              ]                                                              }                                +                                    Nr              k                        .                                              (A10)            
Expression (A10) may be used as the basis for a shift estimation method. Let g be a signal whose shift relative to template signal f is unknown. Let ĝ denote the discrete Fourier transform of signal g. Substituting the transform ĝ for the transform {circumflex over (f)}d in expression (A10) gives                     d        =                                            N                              2                ⁢                                  xe2x80x83                                ⁢                π                ⁢                                  xe2x80x83                                ⁢                k                                      ⁢                          {                                                Arg                  ⁡                                      [                                                                  g                        ^                                            ⁡                                              (                        k                        )                                                              ]                                                  -                                  Arg                  ⁡                                      [                                                                  f                        ^                                            ⁡                                              (                        k                        )                                                              ]                                                              }                                +                                    Nr              k                        .                                              (A11)            
Since the integer value r is not known, there may be ambiguity in the value of shift parameter d. For example, the quantity Nr/k, and therefore, the right hand side of expression (A11), may attain k distinct values modulo N depending on the value of integer r. However, if k equals one, the shift estimation value d is determined unambiguously. Thus, another method for estimating the shift between signals f and g may be summarized by the formula                     d        =                              N                          2              ⁢                              xe2x80x83                            ⁢              π                                ⁢                      {                                          Arg                ⁡                                  [                                                            g                      ^                                        ⁡                                          (                      1                      )                                                        ]                                            -                              Arg                ⁡                                  [                                                            f                      ^                                        ⁡                                          (                      1                      )                                                        ]                                                      }                    ⁢                                    (                              mod                ⁢                                  xe2x80x83                                ⁢                N                            )                        .                                              (A11)            
It is noted that this phase-difference method is much more efficient computationally than the exhaustive cross-correlation and FFT-based cross-correlation methods described above as shown in FIG. 1. The computational effort for the phase-difference method is order N since only one DFT coefficient of the signal g needs to be computed. However, the performance (estimation accuracy) of the phase-difference method degrades at much lower noise levels than the cross-correlation methods as shown in FIG. 2.
Furthermore, the phase-difference method and the cross-correlation methods provide no mechanism for improving estimation accuracy based on knowledge of the noise process which operates on the shifted template f to generate signal g.
Thus, there exists a substantial need for a system and method which could estimate the shift between two signals with a combination of computational efficiency and noise immunity. Furthermore, there exists a need for a shift estimation system and method which could take advantage of information about the underlying noise process to improve estimation accuracy.
There is an implied tradeoff between computational efficiency, and noise resistance as shown by FIG. 2. Thus, there exists a need for a method which could explicitly determine the tradeoff.
The present invention comprises a system and method for estimating the shift between two signals. The shift between the signals may be interpreted as time-delay, spatial shift, rotation angle, etc. For example, the present invention may be used to estimate the shift in a signal g relative to a template signal f. The shift estimation system method comprises: (a) receiving a first signal g, where the first signal g may be represented as a vector g having N components; (b) projecting the vector g to a space with dimension K less than N to obtain a projection vector X having K components; (c) computing measures of distance between the projection vector X and vectors in a set of stored vectors; (d) determining a stored vector p in the set of stored vectors with a minimum distance to the projection vector X based on the measures of distance. The stored vectors are generated from a template signal f, also represented as a vector with N components, by projecting shifted versions of the template signal f to the space of dimension K. The shifted versions of the template signal f may be referred to as shifted template vectors, or simply, shift vectors. The shift estimation method may provide a shift value based on the shifted template vector which generates the stored vector p, where the shift value is as an estimate for the shift between the received signal and the template signal f. The shift value defines the amount by which the template signal f must be shifted to obtain the shifted template vector.
In the preferred embodiment of the shift estimation method, the shifted versions of the template signal f are projected to the space of dimension K by matrix multiplication. In other words, each of the shifted templates are multiplied by an appropriately chosen matrix C. In addition, the projection of the vector g is achieved in the same way, i.e. by multiplying the vector g with the matrix C.
The matrix C may be computed or estimated in a variety of ways. In the preferred embodiment, the matrix C is computed by performing an iterative search in the set of matrices with N rows and K columns in order to maximize an objective function. A simulated annealing algorithm may be used to perform the iterative search.
In one embodiment, the objective function may be evaluated on a candidate matrix Y by (1) multiplying the set of shifted template vectors by the candidate matrix Y to generate a set of image vectors; (2) computing a minimal distance between the image vectors; (3) computing a noise object size associated with the candidate matrix Y; (4) and computing a ratio comprising the minimal distance and the noise object size.
In addition, the method of the present invention may advantageously provide increased estimation accuracy by computing the matrix C based on knowledge of the underlying noise process. For example, given a noise model g=U{f,d} which describes the random generation of signal g in terms of template function f and shift parameter value d, the present invention contemplates (i) multiplying the set of shifted template vectors by a candidate matrix Y to generate a set of image vectors; (ii) computing, for each of the image vectors, a distance to a nearest neighbor in the set of image vectors; (iii) estimating a noise radius (or noise object size) for the distribution of projection vector gY based on the noise model for each value of the shift parameter d; (iv) computing, for each of the image vectors, the ratio of the distance to the nearest neighbor and the corresponding noise radius, (v) determining the maximum of the ratios over the set of image vectors, (vi) iteratively performing (i) through (v) for a series of matrices Y to locate the matrix Y=C which minimizes the maximum of the ratios.
Alternatively, the matrix C may be computed as follows. First, compute the discrete Fourier transform of the template signal f. Second, identify a collection of K/2 distinct integers in the range from 1 to Nxe2x88x921 which maximize the sum of magnitudes of the Fourier transform evaluated on the collection of K/2 distinct integers subject to a spanning condition. Finally, for each integer m in the collection of K/2 distinct integers, assign the cosine and sine waveforms with frequency m to a corresponding pair of columns of matrix C. In one embodiment, the spanning condition is that the greatest common divisor of N and the collection of K/2 distinct integers is one. In another embodiment, the spanning condition is that at least one of the K/2 distinct integer is relatively prime to N.
Let fd denote the vector obtained by circularly shifting the entries of template signal f through d positions. Because template signal f has N components, there are exactly N distinct circularly-shifted versions of template signal f denoted f0, f1, f2, . . . , fNxe2x88x921. The set of stored vectors mentioned in (c) may be generated by projecting any subset S of the circularly-shifted versions f0, f1, f2, . . . , fNxe2x88x921 into the K-dimensional space. The subset S may comprise all of the circularly-shifted versions f0, f1, f2, . . . , fNxe2x88x921.
However, if it is known that the signals g are never generated with values of the shift parameter d in a particular range, the shifted versions fd corresponding to this particular range may be excluded from the subset S. The computational effort for the distance computations (or distance measurement computations), i.e. step (c) above, is proportional to the size of the subset S. Thus, any reduction in the size of subset S advantageously increases the computational efficiency of the shift estimation method of the present invention. Furthermore, it is observed that the average distance between the stored vectors increases as the number shifted versions fd in subset S decreases. This implies a higher probability of correct shift estimation. Thus, knowledge of the type or kind of shift, or the range of possible shifts, may be used to provide increased efficiency and estimation accuracy.
The shift estimation method may be implemented in software configured for execution on one or more processors. Alternatively, the shift estimation method, or any portion thereof, may be implemented by dedicated hardware such as a Field Programmable Gate Array (FPGA).
In the preferred embodiment, the shift estimation system operates on a series of signals which are shifted and perhaps noise-perturbed versions of the template signal f. Each of the signals may be presented as a vector g with N components. The shift estimation method estimates the shift of each of the signals with respect to the template signal f.
The present invention also contemplates a signal analyzer system which comprises a processor coupled to a memory. The processor is configured to execute a shift estimation routine stored in the memory. The shift estimation routine may be configured to implement the shift estimation method of the present invention.